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Summary 


This  project  initiated  various  aspects  of  an  ongoing  study  of 
numerical/analytic  techniques  for  the  identification  of  periodic 
solutions  to  functional  differential  equations .  The  techniques 
developed  apply  to  very  general  classes  of  equations,  and  have  been 
implemented  on  a  variety  of  specific  model  problems. 

’‘Local  *  techniques  refer  to  methods  that  apply  to  the  problem  of 
analyzing  the  Hopf  bifurcation  structure  of  small  periodic  orbits 
of  multiparameter  systems.  A  FORTRAN  code,  BIFDE,  was  written  to 
analyse  generic  bifurcations  of  general  systems  with  infinite 
delay. 

"Global”  tracking  methods  have  been  developed  to  study  the  growth 
and  parameter  dependence  of  global  Hopf  bifurcations. 

Investigations  have  centered  on  the  development  of  spline-based 
approximation  techniques  and  their  implementation  in  a  FORTRAN  code 
FDETRAK. 


Research  Objectives 


As  the  beginning  of  a  large  on-going  project  on  numerical 
techniques  for  the  analysis  of  general  parameter-dependent 
functional  differential  systems,  emphasis  was  placed  on  the 
establishment  of  the  numerical  characteristics  of  some  of  the 
proposed  methods.  An  algorithm  of  the  PI  for  the  analysis  of  the 
bifurcation  structure  of  Hopf  bifurcations  chosen  to  be  implemented 
in  a  general  purpose  code  BIFDE.  A  graduate  student,  Archana 
Sathaye  was  supported  by  the  grant  to  assist  in  the  implementation. 

Both  Fourier  series-based  and  spline-based  methods  are  potentially 
valuable  in  the  approximation  of  the  periodic  solutions  (both  large 
and  small)  of  general  functional  differential  systems.  A 
particular  model  of  neuron  firing  was  chosen  to  investigate  the 
feasibility  of  the  former,  while  spline-based  methods  formed  the 
basis  behind  a  more  general  purpose  FORTRAN  code,  FDETRAK. 
Comparison  of  the  two  methods,  as  well  as  benchmark  studies  of  the 
latter,  formed  the  bulk  of  this  aspect  of  the  project.  A  supported 
graduate  student,  Toanhung  Doanvo,  assisted  in  this  portion  of  the 
project . 

Research^Status 

Local  Analysis:  As  suggested  in  the  original  proposal,  a 
special  case  of  the  Lyapunov-Schmidt  based  stability  procedure 
derived  previously  by  the  PI  has  been  successfully  implemented  in 
a  FORTRAN  code  BIFDE,  whose  development,  coding  and  testing 
constitute  the  major  portion  of  the  Master's  Thesis  of  Archana 
Sathaye  (reference  [3]  below),  completed  July,  1986  at  the 
Virginia  Polytechnic  Institute. 

The  package  is  modular,  and  consists  of  several  routines  which 
perform  one  or  more  tasks.  In  conjunction  with  the  routines 
available  in  the  package,  the  user  is  required  to  provide  a  few 
routines  that  describe  the  specific  system  under  study  (eg, 
bifurcation  data,  characteristic  equation,  quadratic  and  cubic 
nonlinearity).  It  should  be  noted  that  the  code  now  allows 
numerical  analysis  of  a  much  wider  class  of  equations  than 
previously  possible .  The  code  has  been  written  for  VAX  computers , 
and  suggests  no  difficulty  in  the  analysis  jf  small  to  moderate 
siced  systems,  given  that  one  can  obtain  tha  prerequisite  data 
called  for  by  the  program. 


The  code  has  been  extensively  tested  on  models  from  the  areas  of 
mathematical  epidemiology  and  genetic  repression.  In  the  first 
case,  this  numerical  confirmed  an  analysis  of  the  PI  that  treated 
a  cyclic  epidemic  model  in  the  setting  of  functional  differential 
equations  with  infinite  delay.  The  code  treats  the  system  as  a 
system  of  differential  equations,  a  more  natural  setting. 

A  class  of  genetic  repression  models  of  Mahaffy  was  considered 
next.  Application  of  the  code  to  this  delay  difference  model  was 


direct,  and  allowed  a  substantial  enlargement  of  the  scope  of  the 
numerical  analysis  over  that  originally  done. 

A  third  application  was  made  to  a  model  of  "chugging"  in  liquid 
propellant  rockets.  An  analysis  of  the  derivation  of  the  linear 
model  of  Tischer  and  Bellman  revealed  nonlinearities  to  be 
attributed  to  viscous  friction  forces  and  fluid  flow  through  the 
fuel  injection  nozzle  (Torricelli’s  law).  These  nonlinearities 
were  retained,  and  the  resulting  model  was  analyzed  with  the  aid  if 
BIFDE.  No  previous  attempts  had  been  made  to  analyse  the  structure 
of  self -oscillations  in  such  a  nonlinear  model.  Our  work  suggested 
the  existence  of  unstable  periodic  orbits,  which  would  be 
consistent  with  experimental  work  cited  in  the  literature. 

Additional  details  are  available  in  the  cited  thesis,  which  also 
contains  a  code  listing.  With  the  completion  of  additional 
testing,  results  of  this  aspect  of  the  project  will  be  prepared  for 
publication. 

Global  Analysis:  Numerical  tracking  methods  have  been  the 
stress  of  a  sizeable  portion  of  the  project.  The  long  term  goal  is 
to  develop  a  portable  automatic  code  for  the  tracking  of  Hopf 
bifurcations  in  single-parameter  systems,  with  identification  of 
stability  and  secondary  bifurcations,  as  well  as  their  evolutions. 

Work  with  A.  Castelfranco  has  been  completed  on  the  application  of 
Fourier  series  approximation  techniques  to  a  model  of  delayed 
feedback  in  neural  systems.  Such  methods,  while  accurate  and 
efficient  in  handling  periodic  solutions  near  point  of  Hopf 
bifurcation,  are  not  particularly  suited  for  accurate  approximation 
of  large  periodic  orbits.  This  observation  motivates  an  alternate 
spline  based  approximation  scheme. 

The  PI  has  developed  a  working  partially  automatic  curve  tracking 
code  FDETRAK  for  the  implementation  of  such  a  scheme.  As  is  the 
case  of  the  code  BIFDE,  on  input  one  must  provide  the  specifics  of 
the  model  under  study  (eg,  bifurcation  data,  linearizations).  The 
code  automatically  selects  stepsize  strategies  to  continue  the 
one-parameter  family  of  orbits.  Floquet  multipliers  are  computed 
by  approximating  the  Poincare  map  associated  with  the  periodic 
orbits.  In  particular,  a  finite  dimensional  approximation  of  the 
phase  space  leads  to  a  finite  dimensional  approximation  of  the 
period  map.  Multipliers  can  then  be  approximated  with  the  aid  of 
standard  eigenvalue  packages  (eg,  IMSL). 

Currently,  the  code  is  designed  to  analyze  a  restricted  class  of 
one  and  two-dimensional  scalar  delay  difference  equations. 
Stabilities  of  periodic  orbits  are  computed,  and  secondary 
bifurcations  from  the  primary  branch  are  identified.  Code 
parameters  set  a  variety  of  algorithmic  variables  such  as  spline 
order,  grid  density,  dimension  of  approximating  phase  space, 
step-size  criteria,  frequency  of  multiplier  calculation,  stopping 
criteria,  etc. 

Initial  work  in  the  Virginia  Tech  VAX  11/785  has  shown  this 


approach  to  be  beyond  the  capabilities  of  such  machines. 
Supercomputer  time  (Cray-2)  was  obtained  from  the  Minnesota 
Supercomputer  Institute,  and  basic  benchmarks  were  performed.  A 
30*.  1  improvement  in  running  time  was  observed  in  comparison  to  the 
same  (unvectorized)  code  on  the  University  of  Minnesota  -  Duluth 
VAX  750.  Such  results  point  to  the  need  for  the  identification  of 
reliable,  yet  less  time-consuming  algorithms,  as  well  as  the 
investigation  of  parallel  algorithms.  Supported  graduate 
student  Tuanhung  Doanvo  (Virginia  Tech)  has  concerned  the  use 
of  subspace  iteration  methods  to  speed  the  multiplier 
calculations.  Work  along  these  lines  continues. 

As  restricted  as  the  code  now  stands,  it  has  been  tested  on  a 
variety  of  delay-difference  models,  and  has  provided  new 
information  about  them.  In  particular,  numerical  work  appears  to 
confirm  a  conjecture  of  Chow  and  Walther  concerning  the  behavior  of 
a  model  from  nonlinear  optics.  A  epidemic  model  of  Mackey  has  been 
considered  and  the  numerical  results  considered  in  comparison  with 
the  work  of  Sternberg  on  the  subject.  Again,  the  algorithm  is 
highly  reliable ,  although  numerically  intensive . 

A  copy  of  the  code  FDETRAK  constitutes  Appendix  I  of  this  report. 


Research  Publications 

1.  Periodic  solutions  in  a  model  of  recurrent  neural  feedback,  by 
A.  Castelfranco  and  H.  Stech,  SIAM  Journal  of  Applied  Mathematics, 
in  press. 

2.  A  numerical  analysis  of  the  structure  of  periodic  orbits  in 
autonomous  functional  differential  equations,  by  H.  Stech,  to 
appear  in  Dynamics  in  Infinite  Dimensional  Systems. 

3.  BIFDE:  A  numerical  software  package  for  the  Hopf  bifurcation 
problem  in  functional  differential  equations,  by  Archana  Sathaye , 
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.360  wO  -  bifurcation  frequency 

1.52194  pO  -  critical  bifurcation  parameter 

40  n  •  number  of  collocation  points 
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1.  Introduction 

Understanding  the  structure  of  periodic  solutions  in  nonlinear, 
autonomous  functional  differential  equations  is  a  problem  that  often 
arises  when  such  equations  are  used  in  the  mathematical  modeling  of 
"real-world"  phenomena.  Knowledge  of  the  existence,  stability,  and 
parameter  dependence  of  such  periodic  solutions  provides  valuable 
insight  into  the  general  dynamics  of  the  system.  Stable  steady  states 
and  periodic  orbits  are  of  particular  interest  since  they  correspond  to 
observable  states  in  the  system  being  nK>deled.  However,  unstable 
steady  states  and  periodic  orbits  are  of  importance  as  well  since 
(through  variation  of  parameters  in  the  model)  these  solutions  can 
themselves  change  stability  and  therefore,  become  "observable". 

Numerical  simulation  of  the  associated  initial  value  problem  often 
provides  evidence  of  the  existence  of  stable  equilibria  and  stable 
periodic  orbits.  However,  it  is  of  limited  value  in  the  study  of  unstable 
solutions. 

Linearization  provides  a  straight-forward  means  of  analyzing 
equilibria  and  their  stability  types.  A  careful  study  of  the  associated 
characteristic  equation  ideally  leads  to  the  identification  of  the  subset 
of  parameter  space  at  which  a  variation  of  the  system  parameters  can 


/ 


induce  a  qualitative  change  in  the  nature  of  solutions  near  the 
equilibriuno.  Generically,  system  parameters  corresponding  to  the 
existence  of  the  characteristic  value  X««0  correspond  to  branch  points 
of  equilibria,  while  the  existence  of  a  complex  conjugate  characteristic 
root  pairs  X  »±iw  correspond  to  the  existence  of  small-amplitude 
periodic  orbits. 

At  parameter  values  of  this  last  type  (Hopf  bifurcation)  there  is 
now  a  straight-forward  technique  for  the  determination  of  the  stability 
and  parameter  dependence  (i.e.,  direction  of  bifurcation)  of  such  orbits 
|7].  Fixing  all  but  one  system  parameter,  it  is  natural  to  ask  how 
variation  of  the  remaining  parameter  effects  the  periodic  orbit,  inducing 
changes  of  stability  and  secondary  bifurcations.  Towards  this  end,  the 
theory  of  global  Hopf  bifurcation  is  valuable  in  identifying  the  available 
alternatives  [2|. 

This  paper  concerns  the  use  of  numerical  methods  (other  than 
simple  simulation  studies)  to  aid  in  the  analysis  of  both  the  local  and 
global  natures  of  periodic  solutions  to  parameter -dependent  autonomous 
periodic  orbits.  Section  2  discusses  a  numerical  implementation  ot  the 
Hopf  bifurcation  algorithm  of  [7|.  Section  3  outlines  the  use  of 
numerical  tracking  techniques  to  determine  certain  information  concerning 
the  global  bifurcation  picture  in  one- parameter  problems.  The 
usefulness  of  such  techniques  is  illustrated  in  Section  3.  where  the 
result  of  the  analysis  of  a  model  of  nerve  firing  are  described. 


2.  Local  Analysis 

Consider  the  differential  equation 

x'(t)-f(u;xj  |2.ll 


in  which  it  is  assumed  that  x«0  is  an  equilibrium  for  all  values  of  the 


I- 

system  parameters  o€R  .  Somewhat  arbitrarily,  we  have  chosen 
f  :  R^xC  R*',  where  C««C((-r,0),  R**)  is  the  usual  Banach  space  of 

continuous  R*^- valued  functions  on  [-r,0];  other  phases  spaces  can  be 
used  as  well.  Given  adequate  smoothness  (which  we,  henceforth,  assume 
without  mention),  we  may  expand  the  right  hand  side  in  series  form 

x'  (t)  -  L(o)x^  +  H2{a;x^,x^)  +  H3(a;x^,x^,x J  •  •  •  ,  [2.2] 

where  L(a)  is  bounded  and  linear  on  C,  and  and 

H3(o;  are,  respectively,  bounded  symmetric  bilinear  and 

trilinear  forms  on  C. 

The  linearized  problem  reads 

y'  (t)  -  L(a)y^,  (2.3l 

which  has  exponential  solutions  y(t)  *  ^  e^^  if  and  only  if 

ixi  -  L(a)  j  *  ^a;X)f  -  0.  l2.4l 

We  assume  the  existence  of  a  critical  parameter  af*or^  at  which  (2.4| 
has  a  nontrivial  solution  (i.e.,  det^a^;X)«»0)  with  X-±ioj  a  purely 

imaginary  root  pair.  Assuming  simplicit  of  the  root  icj,  it  is  known 
that  the  corresponding  characteristic  vector  ^  is  uniquely  defined  up 
to  a  scalar  multiple.  Furthermore,  the  implicit  function  theorem  shows 
that  there  exists  a  unique  smooth  family  X(o)  of  characteristic  roots 

defined  in  a  neighborhood  of  in  R  and  satisfying  X(a^)«»iu;. 
For  a  near  o^,  we  write  X(a)*p(a)4>iu;(a)  and  For 

simplicity,  we  assume  that  at  there  are  no  other  purely  iautgin^ry 
root  pairs. 

Define  f  "f  (ci)**0  to  be  any  solution  of 


For  X  near  X(a),  let 


« 

^  (a)^o;X(a))«>0  for  a  near  a^. 

A  ^  «  « 

f(o;X)s^  /[f  A' (a;X)f ],  where  A'  — dA/5X.  By  simplicity  of 

the  characteristic  value  X(a)  the  denominator  above  is  nonzero. 

The  following  theorem,  whose  proof  may  be  found  in  [7],  reduces 
the  problem  of  analyzing  the  existence  of  small  periodic  solutions  with 
frequence  near  w  to  that  of  considering  a  scalar  "bifurcation  function". 

Theorem  2.1:  Under  the  above  hypotheses,  there  are  smooth  functions 

G(a;c,i/)  (C- valued)  and  x(t,a;c,i/)  (R^'-valued  and  2  jt/i/- periodic 

in  t)  defined  in  a  neighborhood  of  (a^,0,u;)  in  R  xRxR  such  that  [2.1] 

has  a  small  2 tt/i/- periodic  solution  x(t)  with  (a,i/)  near  (or^.w)  if 

and  only  if  x(t)»x(t,a;c,i/)  up  to  phase  shift,  and  (a,c,i/)  solves  the 
bifurcation  equation 

G(a;  c,  u)  -  0.  [2.5] 

Moreover, 

x(t,Qf;c,i/)  -  2  Re{^ (a)e‘*^**’}c  +  0(c“),  (2,6) 

G  is  odd  in  c,  and  has  the  expansion 

G(a;  c.  i/)  -  (X  -  ui)c  +  M3(a;  u,  X)c^  +  O(c^),  [2.7] 

where  X-X(a), 

Njlor;  v)  m  ^  +  2H2(^„\2  2*"*")  A,  q), 

with  ^(s)— f  (a)e'*^*  for  sSO  and  A.,  A,  q  the  unique  solutions  of 
^Qt;  2i/i)A<,  ,  -  Ho(s»'). 


A(a:0)A2o 

respectively. 


12.8) 


The  imaginary  part  of  [2.6j  can  be  easily  solved  (e.g.»  by 

o 

iteration)  to  obtain  t/<-u;(a)-*-0(c*).  Upon  substituting  this  into  the  real 
part  of  [2.5|,  one  obtains  the  "reduced"  bifurcation  equation 

0  -  g(a;  c)  s  /i(c»)c  +  KjCojc®  +  O(c^).  (2.9l 

Given  K^(a^)#0,  (the  so-called  "generic"  case),  one  can  show  the 

existence  of  nonzero  solutions  c»c  (a)  for  values  of  a  near  for 

which  sgn{/i(a)}«-sgn{K^(a^)}.  If  in  the  case  k»l  /i(a)  increases 

with  a  and  K^(a^)<0,  the  solution  of  [2.9]  near  c— 0  requires 

^(a)>0  (supercritical  bifurcation).  Similarly  K2(a^)>0  corresponds 
to  subcritical  bifurcation. 

Concerning  the  stability  of  the  associated  periodic  orbits,  it  is 

known  [7|  that  if  all  other  characteristic  roots  have  negative  real  parts, 

* 

then  the  periodic  orbits  posses  the  same  stability  type  as  that  of  c 
when  viewed  as  an  equilibrium  solution  of  the  scalar  ordinary 

differential  equation 

c'  -  g(o;  c)  |2.10) 

Thus,  <  0  corresponds  to  an  orbitally  asymptotically  stable 

periodic  orbit,  while  K^(a^)>0  corresponds  to  an  unstable  periodic  orbit. 

The  above  algorithm,  although  usually  too  involved  to  allow  an 
algebraic  determination  of  the  structure  of  Hopf  bifurcations,  does  lend 
itself  to  numerical  evaluation.  This  has  been  recently  implemented  in 
the  FORTRAN  code  BIFDE  by  A.  Sathaye  |6j.  It  is  there  presumed 


that  one  can  obtain  from  the  linearised  equation  |2.3|  analytic 
expressions  for  ^a;X),  as  well  as  the  partial  derivatives  of  ^a;X) 
with  respect  to  X  and  some  (use-chosen)  coordinate  of  a.  It  is  also 
assumed  that  one  is  able  to  identify  critical  values  of  the  parameter 

of^  for  which  a  simple  purely  imaginary  root  pair  X  «±(ji  exists, 

and  the  other  spectral  assumptions  listed  above  hold.  Finally,  it  is 
expected  that  ^^d  can  be  evaluated, 

where  each  of  the  arguments  are  of  the  form  ^j(s)— we**  for 

complex  values  z  and  complex  n- vectors  w. 

Given  the  above  data.  BIFDE  coordinates  the  calculation  of  the 

• 

left  and  right  characteristic  vectors  ^  and  ^  (by  inverse  iteration), 
identification  and  solution  of  the  linear  systems  [2.8]  (by  Gauss 
elimination  with  implicit  pivoting)  and  the  evaluation  of  (hence, 

and  Kj).  The  program  uses  the  partial  derivatives  of  A(a;  X)  to 

compute  (a^)  the  partial  derivative  of  fi  with  respect  to  a 

user -chosen  coordinate  of  o.  and  thereby  determine  the  direction  of 
bifurcation  vnth  respect  to  that  coordinate  of  o. 

The  program  is  complementary  to  BIFDD  of  Hassard  (5|  in  that 
BIFDD  assumes  [2.1]  to  be  of  delay-difference  form,  yet  identifies  the 
required  higher  order  terms  numerically.  In  [5|,  rather  than  making  use 
of  Theorem  2.1,  the  stability  and  direction  of  bifurcation  is  determined 

t 

by  center  manifold  approximation  techniques  and  the  Poincare  normal 
form. 

Remark  2.2:  Given  BIFDE  or  BIFDD.  a  principle  difficulty  lies  in  the 
determination  of  the  bifurcation  data.  That  is,  the  critical  va]ue(s)  of 


/ 


the  system  parameters  o  and  the  associated  frequency  w  of 
bifurcating  periodic  orbits.  For  one- parameter  problems  (k*l),  solution 
of 

del  ^a;  iu;)  —  0  (2.11] 

can  be  obtained  by  standard  root  finding  techniques  (e.g.,  Newton  on 
Quasi-Newton  methods)  provided  the  size  of  the  system  n  is  not 
prohibitively  large,  and  a  sufficiently  accurate  approximation  to  the 
bifurcation  data  is  known  in  advance. 

For  k  large,  one  can  seek  bifurcation  data  by  considering  the 
associated  nonlinear  minimization  problem: 

Minimize  |  A(a;iw)  |  subject  to  the  constraint  that  a  and  w 

lie  within  a  compact  interval  of  u;  values  and  a  lies  within  a 
compact  subset  of  admissable  system  parameters. 

Precise  approximations  for  this  minimization  problem,  although  useful  in 
specifying  the  above  constraints,  are  not  necessary. 

For  k-*2,  one  expects  the  underdetermined  system  [2.1 1]  to  have 
a  one-parameter  family  of  bifurcation  data.  Given  one  set  of 
bifurcation  data  (perhaps  by  considering  the  above  minimization  problem) 
one  can  apply  now  standard  continuation  techniques  to  identify  curves 

2 

in  parameter  space  (a  subset  of  R  )  along  which  {2.11]  has  a  solution. 
Indeed,  in  many  instances,  this  family  of  critical  values  can  be 
parameterized  in  terms  of  ^  itself.  A  simple  example  serves  to 
illustrate  the  point. 

Example  2.3;  Consider  an  ordinary  differential  equation 

x' (l)  -  f(x(i));  X  €  R“  [2.121 

in  which  one  coordinate  Xj  of  x  is  thought  to  act  as  a  feedback  in  one 
of  the  n  equations  in  |2. 12|.  In  studying  the  effects  of  time  delay  in 


this  feedback,  one  replaces  the  appropriate  term  Xj(t)  vdth  Xj(t-r).  To 

determine  the  stabilizing/destabilising  effect  on  equilibria,  one  encounters 
a  characteristic  equation  of  the  form 

p(X)  +  q(X)se~''^  -  0,  (2.131 

where  p  and  q  are  polynomials,  r  corresponds  to  the  length  of  the  time 
delay,  and  s  represents  a  measure  of  the  strength  and  type  (positive  or 
negative)  of  the  feedback.  One  can  algebraically  solve  for  r  (then  s)  in 
terms  of  X  by  considering  [2.13]  and  its  conjugate.  The  details  are 
elementary  and  omitted.  Observe  that  this  provides  a  convenient 
reparameterization  of  [2.1]  in  terms  of  X  rather  than  r  and  s  in  which 
(generically)  with  Re{X}— 0,  Im{X}  determines  the  location  of  on  the 

o 

imaginary  root  curves  in  R",  and  with  lm{X}  fixed,  Re{X}  determines 
the  stability  of  the  equilibrium. 


3.  Global  Analysis 

Consider  [2.1]  in  the  special  case  when  k«l,  and  suppose  that  at 
some  critical  value  of  the  parameter  or^  the  equation  has  been  shown  to 

satisfy  the  hypotheses  of  Theorem  2.1.  Equation  [2.6]  provides  an 
asymptotic  estimate  of  the  resulting  one^parameter  family  of  periodic 
orbits  bifurcating  from  the  equilibrium.  We  discuss  in  this  section 
numerical  methods  for  continuation  of  this  one-parameter  family  away 
from  the  equilibrium,  calculation  of  the  stability  of  the  orbits,  and 
identification  of  secondary  bifurcation  points. 

The  numerical  approximation  of  periodic  orbits  must  not  rely  on 
the  stability  type  of  the  orbits  if  a  complete  global  bifurcation  picture 
is  to  be  obtained.  For  that  reason,  periodic  solutions  are  viewed  as 
solutions  of  a  boundary  value  problem  of  the  form 


13.1] 


F(of,  T,  x^)  -  0 

x(t+l)^x(t),  where  F  :  RxR'^xCxC  R*'.  The  independent  variable  t 
has  been  scaled  so  that  T* periodic  solutions  of  (2.1j  correspond  to 
1 -periodic  solutions  of  |3-l]. 

For  periodic  solutions  x(t)  of  [3.1]  we  introduce  the  finite 
dimensional  approximation 

-  Z  C.  13.2] 

j-1  J  J 

where  the  0 j  represent  appropriate  scalar  1  -  periodic  basis  functions 

and  the  c .  are  in  R*'.  Both  truncated  Fourier  series  and  order 

J 

periodic  B-splines  are  examples  of  approximations  of  this  type.  See  (ij 
and  [3]. 

Collocation  provides  one  means  of  computing  the  coefficients  Cj. 
That  is,  for  N  distinct  nodes  tj  chosen  from  [0.,1.)  one  considers  the 
nN  equations 

F(a.  T.  x(N)(i.))  *  0  [3.31 

in  the  nN+2  unknowns  c^;  o,  and  T. 

VVe  adjoin  to  these  equations  a  scalar  phase  constraint  to  remove 
the  indeterminacy  due  to  the  fact  that  the  phase  shift  of  any  periodic 
solution  of  [3.1{  is  also  a  periodic  solution.  In  the  case  of  trucated 
Fourier  series,  this  corresponds  to  simply  setting  one  of  the  coordinates 
of  the  primary  Fourier  coefficients  equal  to  zero.  However,  there  are 
more  sophisticated  methods  of  instituting  such  a  constraint. 

Finally,  we  adjoin  a  scalar  equation  in  order  to  (in  a  sense) 
specify  which  of  the  one- parameter  family  of  periodic  solutions  is  to  be 


computed.  More  precisely,  given  that  a  solution  to  (3.3| 

has  been  obtained,  one  seeks  a  solution  ^ 

given  arclength  away. 

Having  obtained  two  points  on  the  one- parameter  family  of 
periodic  orbits,  one  can  linearly  extrapolate  (^predict")  an  initial 
approximation  to  the  next  desired  member  of  the  family,  then 
iteratively  improve  that  approximation  ("correct")  by  solving  the  above 
nN+2  simultaneous  nonlinear  equations  by  some  Newton-like  scheme. 

It  should  be  remarked  that  in  the  case  of  ordinary  differential 
equations,  the  use  of  B-splines  has  an  advantage  over  truncated  Fourier 
series  in  that  the  Jacobian  matrices  encountered  are  sparse;  the  precise 
structure  being  dependent  only  on  the  order  of  the  splines  in  use.  For 
functional  differential  equations  such  sparsity  is  lost,  with  the  structure 
of  the  Jacobian  dependent  on  the  form  of  the  equation  [3.1]  as  well  as 
the  parameters  T  and  a.  Despite  this  fact,  splines  possess  certain 
numerical  characteristics  which  speak  in  favor  of  their  use  over 
truncated  Fourier  series. 

Having  computed  an  approximation  x^^^  to  [3.1]  at  the  parameter 
values  a  and  T,  one  determines  the  stability  of  the  orbit  (and 
identifies  secondary  bifurcation  points)  by  computing  approximations  to 

the  orbit’s  Floquet  multipliers.  For  equations  with  finite  delay,  some 

/ 

iterate  of  the  (linearized)  Poincare  map  is  compact  [4|.  The  Floquet 
multipliers  are.  therefore  eigenvalues  of  finite  multiplicity  with  zero  as 
their  only  cluster  point. 

Let  '  denote  an  M  dimensional  approximation  to  the  phase 
space  X  in  use.  VVe  assume  X^^^^cX  and  let  :  X-*X^^^^  denote  a 

projection  of  X  onto  X^^^.  The  approximate  (linearized)  Poincare  map  is 


/ 


defined  to  be 

p(M)  .  p(M)  o  n  I  [3.4l 

where  11  is  the  period  1  map  defined  by  the  linearized  equation 

associated  with  [3.1].  The  eigenvalues  of  serve  as  approximations 

to  the  Floquet  multipliers  of  the  periodic  solution  to  [3.l|. 

Finally,  we  remark  that  due  to  the  autonomous  nature  of  [3.1],  1 
is  always  a  Floquet  multiplier  (4j.  This  fact  provides  a  useful  monitor 

of  the  overall  accuracy  of  the  periodic  solution  approximation  x^^^  and 
the  multiplier  approximation  scheme  described  above. 


4.  An  Example 


We  conclude  with  a  brief  description  of  the  results  of  applying 
the  methodology  described  in  the  previous  sections  to  a  model  in 
physiology.  We  refer  the  reader  to  [1]  for  details,  and  seek  only  to 
indicate  the  kinds  of  information  that  can  be  obtained  when  applying 
these  ideas  to  a  particular  mathematical  model. 

The  two  dimensional  delay-difference  system 
v'  -  h(v)  -  w  +  ,i(v(t  -  r)  -  Vq] 

w'  —  p(v  +  a  -  bw)  [4.1] 

3 

arises  as  a  model  of  recurrent  neural  feedback.  Here,  h(v)— v-v  /3, 
p>0  is  small,  0<b<l,  l-2b/3<a<l,  and  v^  is  the  v  coordinate  of 

the  unique  equilibrium  {vq,Wq)  that  exists  for  [4.l|  when  ^—0,  Fixing 

pt  a  and  b,  one  can  consider  the  associated  Hopf  bifurcation  problem 
in  the  two  remaining  parameters  r  and  /i,  which  are  restricted  to  be 
positive  and  negative,  respectively. 

Linearizing  about  the  equilibrium  of  [4.1],  one  obtains  as 


characteristic  equation 

0  »  X ‘'+(pb-£r)X  +  /7(l  +  b<r)-(/iX+^/?b)e"^ [4.2] 
where  (7=*h'  (v^).  As  indicated  in  Section  2,  one  expects  for  this  two 

parameter  problem  that  there  should  be  curves  of  critical  parameters  at 
which  [4.2]  possesses  purely  imaginary  root  pairs  X«‘±wi.  Since  [4.2| 
has  the  form  [2.13],  one  expects  these  "imaginary  root"  curves  to  be 
parameterized  by  frequency  u. 

Figure  4.1  shows  a  few  of  these  curves  for  p— .08,  a— .7  and 
b=.8.  Along  them  one  is  able  to  determine  the  stability-determining 
constant  K ^  by  numerically  implementing  the  algorithm  discussed  in 

Section  2.  Solid  lines  correspond  to  <  0,  while  dashed  lines 

correspond  to  K ^>0.  One  can  show  that  at  nonintersection  points  of 

the  imaginary  root  curves  the  required  spectral  hypotheses  hold,  and 
that  for  small  /i  that  all  characteristic  roots  must  have  negative  real 
parts.  Thus,  [4.1]  supports  both  stable  and  unstable  periodic  orbits. 

If  one  additionally  fixes  r,  one  can  apply  numerical  tracking 
techniques  similar  to  those  discussed  in  Section  3  to  the  resulting 
one-parameter  problem.  Figure  4.2  depicts  the  global  bifurcation 

diagram  with  r=25.  Solid  lines  indicate  stable  periodic  orbits,  while 
dashed  lines  correspond  to  unstable  periodic  orbits.  See  [l]  for  details. 
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